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Abstract 



o 



In this paper we review air shower data related to the mass composition of cosmic rays above 10 15 eV. After explaining the basic 
relations between air shower observables and the primary mass and energy of cosmic rays, we present different approaches and 
results of composition studies with surface detectors. Furthermore, we discuss measurements of the longitudinal development of 
air showers from non-imaging Cherenkov detectors and fluorescence telescopes. 

The interpretation of these experimental results in terms of primary mass is highly susceptible to the theoretical uncertainties 
of hadronic interactions in air showers. We nevertheless attempt to calculate the logarithmic mass from the data using different 
hadronic interaction models and to study its energy dependence from 10 15 eV to 10 20 eV. 
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1. Introduction 

Knowledge of the mass composition of cosmic rays is of 
key importance for solving the long standing puzzle about the 
origin of high-energy cosmic rays. The mass (and therefore 
charge) distribution can provide strong constraints on the ac- 
celeration of cosmic rays and on their propagation through the 
galactic and extragalactic Universe. Of particular interest are 
measurements of the mass composition in vicinity of structures 
observed in the energy spectrum of cosmic rays, most promi- 
nently the "knee" and the "ankle" at approximately 10 155 eV 
and 10 185 eV, respectively, as well as at energies above 10 195 eV 
at which a flux suppression — possibly the GZK-effect — has 
been observed. Correlated changes of the energy spectrum and 
mass composition can provide important clues to the origin of 
these features. For instance, a rigidity dependent cut-off [1] of 
the spectra of cosmic rays with different charge can lead to a 
gradual increase in the average mass of galactic (see [2] and 
references therein) and extra-galactic cosmic rays [3-5]. On 
the other hand, the interpretation of the ankle and flux suppres- 
sion as a signature of propagation effects [6, 7] of ultra-high 
energy extragalactic protons would require a very light compo- 
sition above the ankle [8]. Furthermore, the location and nature 
of the transition from galactic to extragalactic cosmic rays (or 
lack thereof [9]) leads to distinct predictions of the energy evo- 
lution of the mass composition of cosmic rays in the energy 
region between 10 17 eV and 10 19 eV (see e.g. [10-15]). 

Up to energies of some 10 14 eV, the cosmic ray composition 
can be measured directly and with only minor, or none dis- 
turbing effects from cosmic ray interactions in the atmosphere, 
by employing stratospheric balloon- or satellite-borne experi- 
ments, respectively. Cosmic rays at these energies are covered 
elsewhere, see e.g. [16] or [17] and [18] in this volume. At 
higher energies, measurements of cosmic rays are only possi- 
ble via observations of extensive air showers (EAS). Amongst 



all EAS measurements, those aiming at reconstructing the cos- 
mic ray mass composition are the most difficult ones. This is, 
because the mass of the primary particle can only be inferred 
from detailed comparisons of experimental observables with 
air shower simulations, with the latter being subject to uncer- 
tainties of hadronic interactions at the highest energies. Fortu- 
nately, the data from the LHC at CERN will allow for detailed 
tests of interaction models up to center of mass energies being 
equivalent to cosmic ray energies of 10 17 eV, but an apprecia- 
ble amount of extrapolation of the properties of hadronic in- 
teractions will still be needed to interpret the air shower data 
up to energies of 10 2() eV. In fact, these model uncertainties 
become a dominant source of systematics of estimates of the 
cosmic ray mass composition at the highest energies. There- 
fore, identifying different EAS observables with sensitivity to 
the primary mass is of great importance as these different ob- 
servables should lead to consistent conclusions about the prop- 
erties of the primary particles. Inconsistent conclusions, on the 
contrary, would be indicative for deficiencies of the employed 
hadronic interaction models. 

High energy photons and neutrinos can provide an indi- 
rect approach to the composition of cosmic rays at the high- 
est energies. This is because photo-pion production via the 
Greisen-Zatsepin-Kuzmin effect [6, 7] leads to emission of 
neutrinos [19] and photons from pion decays in the reactions 
P+Jcmb -> P+n° -> P+TY> or likewise p + y C MV -=>n + n + -> 
n + n + + Vfj — > n + e + + v e + + v^. The discovery of cosmo- 
genic photons and neutrinos would prove the existence of the 
GZK energy loss mechanism independently from the cut-off in 
the cosmic ray energy spectrum and would be the unambiguous 
consequence of a light cosmic ray composition at ultra-high en- 
ergies. In the case of heavy primaries, the photo-pion produc- 
tion sets in at much higher energies, because the energy thresh- 
old for this reaction depends on the energy per nucleon and not 
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the total energy. Instead, photo-disintegration is the dominant 
process [20] for energy losses, giving rise to neutrinos from 
neutron decays, but yielding much lower neutrino fluxes in the 
EeV range compared to the photo-pion production by primary 
protons. 

This paper is organized as follows: In the next section we will 
provide an overview of the most important air shower observ- 
ables with respect to measurements of the nuclear composition. 
In Sec. 3 experimental results covering the energy range from 
the knee to ultra-high energies will be presented and searches 
for high energy photons and neutrinos will be briefly discussed 
in Sec. 4. In this review we focus on measurements of the cos- 
mic ray composition from air shower experiments only. Other 
aspects of cosmic rays above the knee are discussed elsewhere 
in this volume [21-23] and in [24—27]. 

2. Air Shower Observables Sensitive to Composition 

There exist a plethora of experimental techniques to charac- 
terize air shower features. In the context of mass composition at 
least two orthogonal measurements are needed to estimate both, 
the energy and mass, of the primary cosmic ray that initiated the 
EAS. This is usually achieved by observing either the longitu- 
dinal development of a shower or by the simultaneous determi- 
nation of the electromagnetic and muonic component of EAS 
at ground level. Further insights can be gained by the study 
of the lateral distribution of particles at ground level, which is 
related to the longitudinal development stage of the shower at 
observation level. 

The ability to deduce the nuclear composition of cosmic rays 
on a statistical basis from these measurements relies to a large 
extent on the theoretical understanding of the shower develop- 
ment and the hadronic interactions that occur within the cas- 
cade. However, the major differences between air showers in- 
duced by different primary masses have a rather simple cause. 
The discriminating power originates from the fact that a pri- 
mary nucleus of mass A and energy E can in good approxi- 
mation be treated as a superposition of A nucleons of energy 
E' = E/A. This superposition model seems plausible because 
the binding energy of nucleons is much smaller than the energy 
of the primary particle. However, the nucleons of a primary 
nucleus are obviously not independent. For instance, the first 
interaction of a nucleus of mass A will happen on average at a 
depth (Xflrst) equal to the interaction length A A < A p and the av- 
erage number of nucleons participating in the first interaction, 
fluctuates. Interestingly enough, it can be shown [28] that 
for (ha) = A Aa/A p (the average number of participating nucle- 
ons from Glauber theory [29]), the inclusive distribution of first 
interactions follows exactly the expectation of the naive super- 
position model, namely that the average depth of interaction of 
the A nucleons is A p . This is the semi-superposition theorem 
that provides a justification of the superposition ansatz under 
more realistic assumptions. 

Whereas the superposition model can explain the main differ- 
ences between air showers induced by different nuclei, it does 
of course not account for nuclear effects such as re-interaction 
in the target nucleus or for the effect of nuclear fragmentation 



on air shower fluctuations (see e.g. [30-32]). More realistic pre- 
dictions of air showers can be obtained by using transport codes 
like Corsika [33], Aires [34] or Cosmos [35] together with 
hadronic interactions models such as Epos [36], QGSJet [37] or 
Sibyll [38] (see [39] for a comprehensive review of air shower 
simulations). 

While these transport codes can give an accurate description 
of air showers given a model for hadronic interactions, it is of- 
ten desirable to qualitatively understand the basic physics be- 
hind a certain air shower observable. For this purpose it can be 
constructive to calculate the air shower development within the 
Heitler-model [40, 41] in which the cascade is approximated 
by a simple deterministic branching model. The original treat- 
ment of electromagnetic cascades has recently been extended 
to hadronic showers in [42-44] and some of these qualitative 
dependencies will be referred to in the following. 

2.1. Longitudinal Development 

The observation of the longitudinal development of the parti- 
cle cascade in the atmosphere is especially well suited for com- 
position studies. In analogy to particle physics detectors, the 
atmosphere acts as a huge homogeneous (as opposed to sam- 
pling) calorimeter in which the air is both the passive material 
that drives the shower development and the active material that 
allows the detection of the produced particles via fluorescence 
or Cherenkov light (cf. Sec. 3.2 and 3.3), as well as the light 
guide through which the signal is transferred to the detector. 

The first inelastic interaction of the primary particle of en- 
ergy E and mass A with an atmospheric nucleus occurs at an 
average depth which is equal to interaction length for inelastic 
nucleus-air collisions, Aa-&x- I n this an d eacn consecutive inter- 
action about one third of the primary energy is transferred from 
the hadronic to the electromagnetic component of the shower 
via decays of neutral pions into photons. The energy transfer 
continues until the charged hadrons decouple from the shower 
by decay into muons and neutrinos, i.e. until their interaction 
length becomes larger than their decay length (see next sec- 
tion). If the charged hadrons start to decay on average after 
«d interaction lengths, then the energy in the electromagnetic 
component is 

-(!)*). (i) 

Since «d > 5 for energies >10 15 eV [45], it follows that most 
of the energy of an air shower can be observed in its electro- 
magnetic part and it is this so-called calorimetric energy which 
allows detectors that can observe the longitudinal air shower de- 
velopment to estimate the primary energy with good accuracy. 

Given this estimator of the energy of the primary particle, 
the orthogonal variable sensitive to its primary mass is the slant 
depth at which the particle cascade reaches its maximum in 
terms of the number of particles, X max . This shower maximum 
is dominated by electromagnetic sub-showers produced in the 
interaction with the largest inelasticity which is usually (though 
not always) the first interaction. An electromagnetic shower of 
energy E reaches its maximum at about 

<^x aS >-^oln£/ £ r (2) 
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Figure 1: Air shower simulation of the shower maximum vs. calorimetric energy. Contour lines illustrate the regions which include 90 % of the showers and the 
inset shows a detailed view at 10 20 eV. 



where Xq ~ 36.62 g/cm 2 is the radiation length in air and 
e£ m « 84MeV is the critical energy in air at which ionization 
and bremsstrahlung energy losses are equal. If the total multi- 
plicity of hadrons produced in the main interaction is N and the 
average hadron energy is E/N, then the shower maximum of a 
primary proton is 
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where both the hadronic interaction length and particle pro- 
duction multiplicity are energy dependent. The factor 2 takes 
into account that neutral pions decay into two photons. Fur- 
thermore, the shower maximum is expected to be influenced 
by the elasticity of the first interaction, K e \. d = E\ ea &IE, where 
£iead is the energy of the highest energy secondary produced 
in the interaction. For interactions with > 0.5 most of the 
primary energy will be transferred deeper into the atmosphere 
and correspondingly the shower maximum will be deeper. We 
are not aware of a consistent treatment of the elasticity within 
a Heitler model for the longitudinal development, however us- 
ing air shower simulations, the dependence on the elasticity fits 
well to 

v , I KelaE 
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The elongation rate [46-48] is a measure of the change of 
the shower maximum per logarithm of energy, 



D = 
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dln£ ' 



For protons and constant elasticity Eq. (4) gives 
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where the changes in multiplicity and interaction length are 
given by 



d\nN 
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Since hadronic interaction models predict an approximately 
logarithmic decrease of A p with energy and N oc E s , D p is 
approximately constant and therefore 



(X p m . dx )*c + D p \nE, 



(8) 



with parameters c and D p being dependent on the charac- 
teristics of hadronic interactions. Using the aforementioned 
(semi-)superposition assumption, one obtains 



(X*) = (X p m . dx (E' = El A)) = c + D p \n(E/A) 



(9) 



and at a given energy the average shower maximum for a mixed 
composition with fractions of nuclei of mass A, is 



(10) 



This equation explicitly demonstrates the relation of (X mdx ) to 
the average logarithmic mass of the cosmic ray composition, 
(lnA> = Y.ifi In A. 

The numerical value of D p from air shower simulations is 
about 25 g/cm 2 (or about 60 g/cm 2 for the change in (X mdx ) per 
decade, D P Q = ln(10) D p ) and therefore proton and iron induced 
air showers are expected to differ by around D p (ln56 - In 1) * 
100 g/cm 2 . Moreover, if the hadronic cross sections and mul- 
tiplicities rise with energy (and if there are no sudden changes 
in the elasticity as for instance suggested in [49]), then Eq. (7) 
leads to Linsley's elongation rate theorem which states that the 



value of D p cannot exceed the radiation length in air, Xq. There- 
fore, Eq. (10) implies that if an experiment measures an elonga- 
tion rate of D > Xq, then a change in the cosmic ray composi- 
tion from light to heavy, d(ln A)/d In E > 0, must be responsible 
for that larger value. 

Results of air shower simulations of (X mllx ) and £c a i are 
shown in Fig. 1. As can be seen, the calorimetric energy is 
indeed a good proxy for the primary energy and exhibits only 
small shower-to-shower fluctuations. And, as expected from 
the relations sketched above, the shower maximum penetrates 
deeper into the atmosphere with the logarithm of the energy and 
is shallower for heavy nuclei than for light ones. The shower-to- 
shower fluctuations in X max are however large and even extreme 
compositions like pure proton vs. pure iron have a considerable 
overlap in their X max -distributions. 

However, these fluctuations carry interesting information 
about the primary particle types and/or the 'mixedness' of the 
cosmic ray composition, that can be experimentally exploited. 
For protons the standard deviation cr of the X max -distribution, 
a 2 = ((^C ax ) 2 ) - (X^ RX ) 2 , is given by the quadratic sum of the 
fluctuation of the first interaction point and the fluctuations of 
the shower development, which in case of the simple Heitler 
model Eq. (4) reduces to 
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Interestingly, for a geometrical treatment of nucleus-nucleus 
scattering [50], the relative fluctuations of the multiplicity, 
<t(N)/N, are constant. Furthermore, if the relative fluctuations 
of the elasticity do not depend on energy, then cr 2 is expected 
to decrease with energy reflecting the logarithmic rise of the 
proton-air cross section. 

For primary nuclei, the naive expectation from the super- 
position model would be cr A = cr p / ^fN. As will be seen in 
Fig. 11, air shower simulations predict a cr p of about 60g/cm 2 
at 10 18 eV and correspondingly <x Fe should be at the level of 
8 g/cm 2 whereas a much larger value of =s 20 g/cm 2 is predicted. 
The reason for this discrepancy is that although according to 
the semi-superposition theorem the positions of elementary 
nucleon-nucleon interactions are distributed oc exp{-X/ A p ), the 
individual positions are not independent. For instance, the n 
nucleons participating in the first interaction have all the same 
X = Xfi rst . These correlations together with the event-by-event 
variance of n are the main reason for cr A > cr p / yN. In ad- 
dition there are two extreme cases for the spectator nucleons: 
Either complete fragmentation (i.e. all A - n spectators propa- 
gate independently after Xf\ lst ) or a single spectator nucleus with 
A' - A - n, where the latter case without fragmentation would 
result to an additional broadening of the X max -distribution. 

In the interesting case of a mixed composition with fractions 
/; of nuclei of mass A,-, the combined X m . ix -distribution follows 
from elementary statistics [51]. If each component has an aver- 
age shower maximum of (X max ),- with width cr, then 

cr(X max ) 2 = (cr 2 ) + (({X mdx } 2 ) - <X max ) 2 ) (12) 

where (X max ) = «X max ),) is the mean of the combined distri- 
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Figure 2: Composition sensitivity of a combined measurement of the shower 
maximum and its fluctuations. Dots denote different models of the extragalactic 
cosmic ray composition at Earth: Pure proton composition in the dip-model [8] 
and mixed composition with a large [11] or small [23] maximum energy of 
the sources and different spectral indices p at the source. Energies range from 
10 18 - 5 eVto 10 20 eV with a spacing of A lg E = 0. 1 . Red lines are simulations 
with Sibyll at the same energies for various two-component transitions. The 
contour defined by these transition contains all other possible mixtures for 1 < 
A < 56. 



bution. In case of a composition comprised of two components 
only, this reduces to 

cr(X max ) 2 = fcr\ + (1 - /) cr 2 + /(l - /)(A<X max » 2 . (13) 

Therefore, depending on the separation of the mean values 
A{X m . dx ) and the fraction / of component 1, it can happen that 
the combined distribution is broader than the individual dis- 
tributions because their separation adds to the total width. In 
terms of average logarithmic mass, Eq. (12) can be rewritten as 

cr{X m , x ) 2 = {cr 2 } + D 2 ((In 2 A) - {In A) 2 ) , (14) 

which demonstrates that <r(X max ) 2 is proportional to the vari- 
ance of In A. Starting from this equation, Linsley proposed to 
analyze air shower data in the cr(ln A)- (In A) plane [52] which, 
however, involves the need to specify the dependence of <x A on 
In A. In this article we will use a similar approach below (cf. 
Sec. 3.3), but stick to the experimental cr(X max ) and (X max ) val- 
ues. For the purpose of comparing data at several energies to 
model predictions, it is useful to subtract the (energy depen- 
dent) values for iron nuclei from both the data and predictions. 
The approximate energy independence of this subtracted vari- 
ables is illustrated in Fig. 2. Two-component transitions as pre- 
dicted from air shower simulations with Sibyll are shown as 
curves that are, as expected from Eq. (13), parabolas in the 
0"(Xnax)-(^max) plane. These parabolas define a closed con- 
tour that contain all possible combinations of mass mixtures 
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Figure 3: Air shower simulation of the number of muons vs. electrons at ground level for a vertical shower observed at 800g/cm 2 . Contour lines illustrate the 
regions which include 90 % of the showers and the inset shows a detailed view at 10 20 eV. 



for A e [1,56]. To illustrate the discriminative power of the 
0"(^max)-(^max) combination, three models for energy evolu- 
tion of the extragalactic cosmic ray composition [8, 11, 23] are 
shown as well. 

2.2. Particles at Ground 

Another way of detecting cosmic rays and to estimate their 
mass is given by the measurement of particle densities of air 
showers at ground. In the calorimeter analogy of the previ- 
ous section, this would correspond to a calorimeter with only 
one active readout plane and correspondingly this measurement 
technique is more susceptible to shower-to-shower fluctuations. 
Nevertheless, ground measurements are still frequently used in 
cosmic ray detectors because of their geometric acceptance and 
high duty cycle. 

An estimate of the qualitative dependencies of the number 
of muons and electrons on primary mass and characteristics of 
hadronic interactions can again be obtained within the Heitler 
model. Given the average multiplicity N of each interaction, 
the energy of charged and neutral pions in a shower initiated by 
a primary proton of energy E is E n - E/N" after the nth in- 
teraction if one (somewhat unrealistically) assumes an energy 
independent multiplicity. This energy splitting continues un- 
til the charged pion energy reaches the decay energy at which 
the hadronic interaction length A[ nt becomes equal to the decay 
length /ld ec - pycT, where p is the height-dependent density of 
air, y denotes the Lorentz -boost and r is the pion lifetime. For a 
shower with incident angle 6 in an isothermal atmosphere with 
scale height ho, the density at slant depth X is 
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(15) 



Therefore, the condition A m - /Idee leads to a decay energy that 
is independent of the interaction length. It is reached after «d 



n d N-" d = 
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and therefore 



«d = -- 



xu I h m n c 2 \nN \ 
yy -l\ CT E COS0J 



laN 



(16) 



(17) 



where W_i denotes the lower branch of the Lambert- W function 
(see e.g. [53]). The decay energy is then given by 



E 



(18) 



for which we find numerical values of a few tens of GeV and a 
slow decrease with primary energy in agreement with the esti- 
mates of [43]. The total number of muons produced in a shower 
is equal to the number of pions with E n 



and therefore 



with 



InfiV 
lnJV ' 



(19) 



(20) 



where the factor | gives the approximate fraction of charged 
pion secondaries. Air shower simulations predict ft to be in the 
range of 0.88 to 0.92 [42], corresponding to effective multiplic- 
ities from 30 to 200 in Eq. (20). It is interesting to note, that 
because the interaction length drops out in the calculation of nj 
(cf. Eq. (16)), the number of muons at ground are expected to 
be independent of /l; nt . 

The number of electrons at shower maximum, i.e. at the point 
at which the electron energies become too low to produce new 
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particles (E e = e™), can be estimated from the total amount of 
energy in the electromagnetic cascade given by the primary en- 
ergy minus the energy in muons. Since 
of electrons is 
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where the last approximation can be made at high energies at 
which the energy fraction transfered to muons becomes small. 

Using again the superposition model and substituting E with 
E' = E/A, one obtains the following relations for nuclear pri- 
maries: 

El A 



N A 



N p 

1 v c,max 



and 



N: 



(22) 



(23) 



So, whereas the number of electrons at shower maximum gives 
a good estimate of the primary energy independent of the com- 
position, the number of muons can be used to infer the mass 
of the primary particle, since it grows with A 1 ^. Moreover, 
the evolution of the muon number with energy, dA^/d In E, is a 
good tracer of changes in the primary composition. Just as in 
the case of the elongation rate of the longitudinal development, 
a constant composition gives dA^/d In E = ft and any departure 
from that behavior can be interpreted as a change of the average 
mass of the primaries. 

Unfortunately, the experimental situation is more compli- 
cated, because surface detectors do not observe the num- 
ber of electrons at shower maximum, but at a fixed depth 
Aground / cos 8. If the detector and shower maximum are sepa- 
rated by AX = X groun( j/ cos - X max , then only the attenuated 
number of electrons is observed with 



N e.ground ~ N e,max exp [ 



AX 



(24) 



where A m 60g/cm 2 is the attenuation length of the number 
of electrons after the shower maximum. Since heavy primaries 
reach their shower maximum at smaller depths than light ones, 
the number of electrons on ground is expected to be compo- 
sition sensitive as well, with a larger electron number for air 
showers initiated by light primaries. This feature is visible in 
Fig. 3, where A^ vs. N e is shown for air shower simulations at 
different energies for a detector located at 800 g/cm 2 . As can be 
seen, the In A^ and lnN e observables are basically rotated from 
the desired quantities, In A and ln£\ Due to the steeply falling 
cosmic ray spectrum, this rotation causes a complication in the 
analysis of air shower data, because showers of equal lnN e 
are enriched in light elements (cf. Sec. 3.1 for a description of 
unfolding methods to overcome this problem). Furthermore, 
Eq. (24) implies that given the X max fluctuations explained in 
the last section, the relative fluctuations of the electron number 
are expected to be quite substantial, 



^"(^f.ground) ^"(^max) 
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These attenuation effects can be reduced considerably by 
choosing an appropriate detector site which is situated at a 
height close to the shower maximum. The exponential atten- 
uation Eq. (24) is only valid far from the maximum, whereas in 
its close vicinity the shower size is nearly invariant under small 
displacements from the maximum (see Fig. 9 below). Since the 
simulations in Fig. 3 were performed at a fixed ground depth of 
800 g/cm 2 , the evolution of the attenuation effect with distance 
to the shower maximum can be seen indirectly: At low energies 
where the observation level is far from the shower maximum, 
the difference in the number of electrons between proton and 
iron primaries is large and diminishes while the shower maxi- 
mum approaches the ground level at higher energies. 

Besides the measurement of the number of electrons and 
muons, experiments with surface detectors have further means 
to determine the shower age (i.e. the distance to the shower 
maximum) by studying the shape of the particle densities with 
respect to the distance to the shower core. These measurements 
of the lateral distribution as well as other additional composi- 
tion sensitive variables from ground detectors will be discussed 
in Sec. 3.1.2. 

2.3. Model Uncertainties 

The physics of air showers is very well understood in terms 
of particle transport through the atmosphere and for electro- 
magnetic showers it is currently believed that they can be 
modeled without any significant uncertainties. In the case of 
hadronic showers, however, there is a fundamental lack of the- 
oretical and experimental knowledge of the characteristics of 
hadronic interactions (see e.g. [55, 56] for recent discussions of 
hadronic interactions in air showers). Since most of the interac- 
tions in an air shower are 'soft interactions', i.e. occur with only 
a small momentum transfer, perturbative QCD is not applicable. 
The phenomenological interaction models that have to be used 
instead, are constrained by low energy experimental data only, 
but even at these energies the full phase space relevant for air 
shower interactions is not fully covered [45, 57-59]. At ultra- 
high energies, the center of mass energies of the first nucleus- 
air interactions are beyond accelerator energies and correspond- 
ingly the models solely rely on extrapolations. 

This somewhat bleak situation is currently alleviated by 
special-purpose programs to measure particle production data 
relevant for cosmic rays (e.g. [60, 61]) and by the wealth of new 
data from the multipurpose detectors at the LHC with which 
the interaction models can be tested and re-tuned. At its maxi- 
mum center of mass energy of 14TeV the LHC will eventually 
reach the equivalent of 10 17 eV in the laboratory system and 
thus cover most of the energy range at which galactic cosmic 
rays are detected with air showers. First comparisons of the 
7 TeV data to interaction models used in air shower simulations 
suggest that the LHC data are bracketed by these models [62]. 

The effect of uncertainties in hadronic interaction charac- 
teristics on air shower observables has been recently studied 
in [54] and [63]. Some of the results from [54] are displayed 
in Fig. 4 where response of air shower observables to ad-hoc 
changes of the Sibyll interaction model are shown. As can be 
seen, all air shower observables are highly susceptible to such 




Figure 4: Sensitivity of air shower observables on characteristics of hadronic interactions as a function of a change in model characteristics growing logarithmically 
from zero at Tevatron energies to the quoted values on the x-axis at 10 19 eV (proton showers at 10 19 5 eV, adapted from [54]). 



changes but to a different degree. For instance, the muon num- 
ber and X lmx are very sensitive to the multiplicity (as expected 
from Eqs. (4) and (20), whereas the X m - dx fluctuations are not. 
And likewise there is a strong dependence of X max and <x(X max ) 
on A but not so for (as suggested by /(-independence of n& in 
Eq.(16)). 

3. Measurements of the Nuclear Composition 

Strictly speaking, no air shower observatory measures the 
primary composition of cosmic rays. Instead, one or more of 
the mass sensitive observables from the last section can be mea- 
sured and the data can then be interpreted in terms of primary 
mass by a comparison to air shower simulations using hadronic 
interaction models. Since different air shower observables react 
differently to changes in the characteristics of hadronic inter- 
actions, one may hope to diminish the model dependence of 
primary mass estimates by comparing the results from different 
observables. 

In this section we will therefore first review measurements 
of composition related air shower observables from particle de- 
tectors at ground and from optical detectors. The data will then 
be compared to each other in terms of the average logarithmic 
mass in Sec. 3.4. 

3.1. Particle Detectors 

3.1.1. N e -N M Method 

Measuring electron and muon numbers (and their fluctua- 
tions) has become the first and most commonly employed tech- 
nique applied to infer the cosmic ray composition from EAS 
data. The basic principles of this approach can be understood 
very intuitively, since the sum of electron and muon numbers 
measured at ground relates to the primary energy, while the ra- 
tio of the muon to electron number relates to the primary mass 



for reasons described in the previous section. These general 
features can also be observed in Fig. 3. Fukui et al. [66] and 
Khristiansen et al. [67] were the first to study the muon num- 
ber fluctuations in the knee region and they were also the first 
concluding an enrichment by heavy nuclei above the knee en- 
ergy [68]. Similar conclusions were drawn in [69] based on a 
larger data set. 

In the most classical approach, electron-muon discrimina- 
tion is achieved by employing a combination of unshielded and 
shielded scintillation detectors at ground level. Recent exam- 
ples include AGASA [70], CASA-MIA [71], EAS-TOP [72], 
GRAPES [73], KASCADE [74], KASCADE-Grande [75], 
Maket-ANI [76], GAMMA [64], and Yakutsk [77]. The Pierre 
Auger Observatory [78] operates water Cherenkov detectors of 
1.2 m depth which also enables limited muon identification. 
The electromagnetic particles (electrons, positrons, and pho- 
tons) are more numerous than the muons but their mean en- 
ergy at detector level is only some lOMeV while that of the 
muons is about 1 GeV. The former thus produce a large num- 
ber of relatively small Cherenkov pulses whereas muons pro- 
duce a small number of large pulses. This type of discrimina- 
tion technique was pioneered at the Haverah Park detector [79]. 
Since FADC readout was not available at that time, signal traces 
of the detectors were photographically recorded from oscillo- 
scopes. IceTop at the South Pole [80] uses tanks of frozen ice 
instead of water for obvious reasons. Sensitivity to the primary 
composition is achieved by measuring the dominant electro- 
magnetic component at ground level in IceTop in coincidence 
with high energy muon bundles in IceCube (muon threshold 
about 500 GeV), originating from the first interactions in the at- 
mosphere. Similar combinations of surface and underground 
detectors have been used e.g. by EAS-TOP and MACRO at the 
Gran Sasso [81] and at the Baksan underground laboratory [82]. 
Such methods are potentially very attractive for composition 
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Figure 5: Unfolded fluxes from GAMMA [64] and KASCADE [65] using two different interaction models. The spectra of elemental groups have been adjusted by 
a common factor (within 35%) to match the all-particle spectra of the two experiments. KASCADE points were slightly displaced in energy for better visibility. 



studies, as the TeV energy muons measured deep underground 
probe the energy per nucleon while the surface array probes the 
energy per particle of the cosmic ray primary. The disadvantage 
is the reduced statistical accuracy of the muon measurements 
and the limited solid angle available for coincident cosmic ray 
observations. Finally, Tibet-ASy [83] complemented their air 
shower array with emulsion chambers and so-called burst detec- 
tors, and Telescope Array (TA) [84] employ solely unshielded 
scintillation detectors and thus do not discriminate muons from 
electrons at detector level. Composition studies in these cases 
rely on additional measurements discussed below. 

The importance and proper treatment of air shower fluc- 
tuations in any analysis of the nuclear mass composition 
has already been emphasized. It is important to realize 
that electron and muon numbers do not fluctuate indepen- 
dently on event-by-event basis, but are mutually correlated (cf. 
Fig. 3). All these effects can be properly accounted for by 
a two-dimensional unfolding technique, first utilized by the 
KASCADE-Collaboration [65, 85]. This approach yields a set 
of energy spectra of primary mass groups, such that their re- 
sulting simulated double differential electron and muon num- 
ber distribution, dPN/dNedN^, resembles the observed one. 
The only input required for solving the underlying mathemati- 
cal equations are the so-called kernel functions describing the 
transformation (E,A) — > {N° h:i ,N° hs ) including their shower- 
by-shower fluctuations. These functions need to be gener- 
ated from air shower simulations (including detector response 
functions) and thus will depend on the chosen hadronic in- 
teraction model employed in the EAS simulations. The level 
of systematic uncertainties imposed by the interaction models 
can be studied by constructing kernel functions using differ- 
ent hadronic interaction models and comparing the correspond- 
ing unfolded results (cf. Fig. 5). Moreover, the observed two- 
dimensional (N° hs ,N° hs ) distribution can be directly compared 
to the forward folded distribution when using the unfolding re- 
sults as input. Based on such kind of systematic studies the 



KASCADE Collaboration concluded that neither QGSJet-01 
nor Sibyll2.1 is able to describe the full range of energies 
consistently. Specifically, a muon deficit (or electron excess) 
has been found for Sibyll2.1 simulations when compared to 
data [65]. Expressed in terms of the mean logarithmic mass 
(see below) the two models differ systematically by about half 
a unit but yield the same basic result of a light (He dominated) 
composition at the knee with a change towards a heavy com- 
position at higher energies. Moreover, the data are consistent 
with the assumption of a rigidity dependent change of the knee 
energy. Similar conclusions about an increasing mass across 
the knee have been drawn e.g. in [81, 86, 87]. More recently, 
the GAMMA Collaboration [64] performed a forward folding 
approach to their electron and muon distributions. In this case, 
specific spectral shapes are to be defined (here: power-law in- 
dices below and above the knee, E^/Z, and the fluxes of four 
different mass groups). By optimizing the spectral parameters 
one tries to match the forward folded mean N e vs. iV„ distri- 
bution to the measured data. The results of this analysis are 
included in Fig. 5 and are shown as full and dashed lines. A 
good description of the all-particle spectrum is again found for 
a light composition at the knee with an increasing fraction of 
heavy primaries above. However, differences between forward 
folding by QGSJet-01 and Sibyll2.1 lead to A(lnA) * 0.7-0.8 
with an overall rather light composition (cf. Fig. 14). Data from 
Tibet ASy [88], on the other hand, seem to suggest an increas- 
ing mass composition already below the knee leading to a heavy 
dominated composition at the knee. This conclusion, however, 
is based on a number of assumptions and a rather limited statis- 
tics of events. Since the set-up and event selection used in 
Ref. [88] was essentially blind to heavy primaries, the heavy 
component has been deduced in Ref. [88, 89] rather indirectly 
by subtracting the flux of the light component (p+He) from the 
all-particle spectrum measured by Tibet-Ill [90]. 

Only very recently, data from the KASCADE-Grande exper- 
iment reached sufficient statistics allowing the unfolding analy- 
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Figure 6: Unfolded fluxes from KASCADE and KASCADE-Grande (open and 
full symbols, respectively) based on the QGSJet-II-02 interaction model (data 
are from [91,92]. 



with the shower age 0.5 < s < 1.5, the Moliere radius r M , and 
the electron shower size N e . A large variety of Lateral Distribu- 
tion Functions (LDFs) is found in the literature (named Greisen, 
Greisen-Linsley, etc.), often being specific modifications of the 
NKG-function aiming at providing an optimum description for 
either electrons or muons, or both, or are simple power-laws be- 
ing optimized to describing densities measured by scintillators 
or water Cherenkov detectors, etc. [96-98]. 

Besides being a tool to reconstruct shower sizes, the actual 
shape of the LDF also contains information about the under- 
lying particle physics in the EAS and, thereby, also about the 
mass of the primary particle. Generally, showers initiated by 
heavy primaries and reaching their maximum at high altitudes 
will exhibit a flatter LDF than those initiated by light primaries 
and developing deeper into the atmosphere. This feature is ob- 
served both for electrons and muons and the steepness of the 
LDF can be extracted for samples of EAS (properly selected 
by shower size and zenith angle) [99] or, if sampling statistics 
does allow, also on event by event basis [97]. A reanalysis of 
data from Volcano Ranch [100] was performed in Ref. [101] 
and yielded an iron fraction of about 75% at E ^ 10 18 eV, if the 
QGSJet-01 model was adopted. The sensitivity of the LDF to 
the primary mass, however, is weaker as compared to the N e -N^ 
method or those techniques discussed below. 



ses to be extended to energies nearing 10 18 eV (cf. Fig. 6). Be- 
cause of poorer detector resolution in KASCADE-Grande, the 
analysis has been limited to only three mass groups. However, 
the preliminary unfolding results [92] confirm earlier findings 
of KASCADE and indicate a very heavy composition at about 
10 17 eV. Similar conclusions about the mass composition us- 
ing electron-poor and electron-rich showers have been drawn 
in [93] with a knee-like structure observed at 8 x 10 16 eV that 
could be identified with the heavy component of primary cos- 
mic rays. 



3.1.2. Steepness of Lateral Distribution Function 

The determination of the lateral density distribution of par- 
ticles in EAS is the most fundamental task of a particle detec- 
tor, because the total number of charged particles, electrons, or 
muons is obtained by fitting the radial fall-off of the measured 
particle densities to proper analytical functions and perform- 
ing the radial integration. The energy and mass of the primary 
particle can then be deduced from the electron and muon num- 
bers, as described in the previous section. Historically, choices 
of parameterizations of both electron and muon lateral distri- 
butions were influenced very much by the seminal review of 
Greisen [94] in which the analytical calculations performed by 
Kamata and Nishimura [95] for electromagnetic showers were 
generalized to the Nishimura-Kamata-Greisen (NKG) function: 



PNKc(r,s,N e ) = 



T(4.5 - s) 



Ne 

r 2 M 27tF(s)Y(4.5-2s) 

■-2 i x s-4.5 

1 + ^1 

rMl \ r M j 



(26) 



3.1.3. Muon Tracking and Timing 

Interest in reconstructing the the mean heights of production 
of muons in EAS to learn about the longitudinal shower devel- 
opment date back to the 60ies of the last century [102]. How- 
ever, at that time neither the angular resolution of EAS arrays 
nor the accuracy of shower simulations were appropriate to the 
use such measurements for investigating the mass composition 
of primary cosmic rays. The technique was revived in the 90ies 
with the Cosmic Ray Tracking (CRT) detectors at HEGRA 
(aiming at tracking of electrons and muons) [103], the Muon 
Tracking Detector at KASCADE [104] or at GRAPES [105]. 
In these detectors the orientation of the muon track is measured 
with respect to the shower axis and the muon production height 
is then reconstructed by means of triangulation. For such de- 
tectors one needs to find a compromise between a sufficiently 
long base-line at the ground (i.e. distance between shower core 
and tracking detector) and a sufficiently large area for the track- 
ing detectors. Increasing the average distance improves the ac- 
curacy of the triangulation but shifts the tracking detector in 
regions of smaller muon densities. A typical compromise at 
energies around and above the knee is about 100 m baseline or 
more and at least 100 m 2 muon tracking area. Based on such 
measurements, no or only a small rise of the mean mass above 
the knee was concluded from CRT at HEGRA [106]. However, 
data were binned as a function of the electron shower size which 
imposes a strong bias towards electron rich showers, i.e. to light 
primaries. Measurements of the muon production height with 
KASCADE-Grande [107] are compatible with a clear transition 
from light to heavy cosmic ray primary particles with increas- 
ing shower energy across the knee. 

An attractive alternative to infer the mean muon production 
depth in EAS was suggested in [108, 109]. Instead of per- 
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forming a geometrical triangulation by complex and expensive 
tracking detectors, one measures the time delay of muons in 
each detector station at ground level with respect to the shower 
front. Since muons are in general produced close to the shower 
axis, those muons produced at high altitude can reach a detec- 
tor station at some distance to the shower core at a shorter av- 
erage path length as compared to muons produced deep in the 
atmosphere. Therefore, the relative time delay to the shower 
front plane is a direct measure of the mean Muon Production 
Depth (MPD) in the atmosphere. Similarly to the aforemen- 
tioned muon tracking detectors, the resolution of measurements 
of the MPD improves both with increasing muon numbers ob- 
served in a detector station and with increasing distance of the 
station to the shower core. However, distortions become large 
close to the shower core because of limited FADC sampling 
times and electron contamination. Thus, a proper minimum 
distance to the shower core needs to be maintained. 

A first application of this novel method to real data has 
been presented by the Pierre Auger Collaboration only very 
recently [110]. From the MPD distribution measured in each 
event, the quantity X'™ x can be extracted by fitting a Gaisser- 
Hillas function [111] to this distribution, where X™ ax can be 
interpreted as the point where the production of muons reaches 
the maximum along the cascade development. The resolution 
in individual events ranges from about 120g/cm 2 at the lower 
energies to less than 50 g/cm 2 at the highest energy and the sys- 
tematic uncertainties on (X™ ax ) were reported 1 1 g/cm 2 . Within 
these uncertainties, a mixed composition with a mild change to 
a heavier composition is reported in [110] for energies above 
2 ■ 10 19 (cf. Fig. 14). 

3.1.4. Rise-Time 

Further information about the primary mass of UHECR is 
encoded in the time profile of shower particles. This was first 
noted by Bassi et al. [112] and was used in the context of pri- 
mary composition at Haverah Park [113-115]. An advantage 
of this technique is that it does not require electron-muon dis- 
crimination at the detector level, but makes use of the fact that 
muons are mostly detected near the shower front and that elec- 
trons are subject to stronger attenuation in the atmosphere than 
muons. Full exploitation of this technique, however, requires 
state-of-the-art FADC read-out electronics. Moreover, it has 
been noted that the signals from water Cherenkov detectors lo- 
cated upstream of an inclined shower exhibit faster rise-times 
than those located downstream. This is called the 'rise-time 
asymmetry' [116]. As rise-time we shall understand here the 
time it takes for the integral signal to increase from 10 to 50% 
of its full value. 

For simple geometrical reasons, no rise-time asymmetry is 
expected for vertical showers. For very inclined showers one 
again finds only very small rise-time asymmetries between in- 
ner and outer detectors. This can be understood because of the 
absence of the electromagnetic component and dominance of 
the high energy muon component in inclined showers, the latter 
being basically free of asymmetries. The zenith angle at which 
the rise-time asymmetry becomes maximum, 6 max , was demon- 
strated to be correlated with the shower development [116] and 



O 

O) 10 2 





-•- Tunka25 event 




photon 


SSSsllllliss 


— proton 








iron 







40 50 60 70 100 200 300 400 500 

core distance [m] 



Figure 7: Example of a lateral Cherenkov light distribution as measured with 
non-imaging Cherenkov detectors. Data points are taken from [118] (E = 
0.12EeV) and compared to ten simulated [33] air showers for three different 
primary particle types using the hadronic interaction model QGSJetII [37]. 

thereby to the average mass of the primary particles. Like the 
Ne-Nfj unfolding, this method cannot be applied on an event-by- 
event basis, which however is no major limitation. A first appli- 
cation to data was presented in Ref. [1 17] indicating a transition 
from a light to heavy composition in the energy range 3 x 10 18 
to4x 10 19 eV (cf. Fig. 14). 

3.2. Non-imaging Cherenkov Detectors 

If atmospheric conditions permit, particle measurements at 
ground are often complemented by observations of Cherenkov 
or fluorescence light from EAS. 

Cherenkov light emitted by extensive air showers in the at- 
mosphere has been known for many years to contain informa- 
tion on shower development [119]. Non-imaging observations 
by operating PMTs with large Winston cones looking upwards 
into the night-sky are perhaps the simplest and most straight- 
forward technique of EAS observations. The basic principles 
have been worked out in [120, 121] and the method was first 
successfully applied at energies around the knee by the AIRO- 
BICC detectors installed at the HEGRA array [122] and by 
CASA-BLANCA [123] in Utah. More recent measurements at 
the Tunka [118, 124] and Yakutsk [125] arrays have extended 
these measurements up to ultra-high energies. 

The lateral distribution of Cherenkov light at ground is the 
result of a convolution of the longitudinal profile of charged 
particles in the shower with the energy threshold for Cherenkov 
emission and the Cherenkov emission angle which both depend 
on the air density and thus height. Moreover, the angular distri- 
bution of electrons in a shower due to multiple scattering con- 
tributes to the lateral extend of Cherenkov light at ground. Be- 
cause the electron energy distribution (and thus the number of 
electrons above Cherenkov threshold) as well as their angular 
distribution are universal in shower age [134-137], the non- 
imaging Cherenkov technique provides a model independent 
method to measure both, the calorimetric energy and shower 
maximum of air showers. 
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Figure 8: Measurements of (X raax ) with non-imaging Cherenkov detectors (Tunka [118, 126], Yakutsk [127, 128], CASA-BLANCA [123]) and fluorescence 
detectors (HiRes/MIA [129], HiRes [130], Auger [131] and TA [132]) compared to air shower simulations [133] using hadronic interaction models [36-38]. HiRes 
and TA data have been corrected for detector effects as indicated by the (A) values (see text). The right panel shows a zoom to the ultra-high energy region. 



A characteristic feature of the lateral light distribution at 
ground is a prominent shoulder at around 120 m from the 
shower core (cf. Fig. 7) which is due to the strongly forward 
beamed emission of the Cherenkov light (6?£ * 1.4°) from near 
the shower maximum in the atmosphere. The slope of the lat- 
eral distribution measured within this 120 m is found to depend 
on the height of the shower maximum and hence on the mass 
of the primary cosmic ray nucleus. The overall Cherenkov in- 
tensity at distances beyond the shoulder, on the other hand, is 
closely related to the calorimetric energy. 

The (X max ) measurements from BLANCA [123], 
Tunka [118, 126] and Yakutsk [125] are shown in Fig. 8. 
At low energies (E < 10 16 eV) the three measurements disagree 
by up to 40g/cm 2 , but all three detectors observed small elon- 
gation rates above 5 x 10 15 eV, indicating a change towards a 
heavier composition. At around 10 17 eV the absolute values of 
(X max ) from Tunka and Yakutsk are approaching the simulation 
results for heavy primaries and beyond that energy the average 
shower maximum increases again towards the air shower 
predictions for light primaries. At even higher energies, only 
the Yakutsk array measured (X max ) with Cherenkov detectors 
and we will discuss this range in the next section together with 
the data from fluorescence telescopes. 

3.3. Fluorescence Telescopes 

After the first prototyping and detection of fluorescence light 
from air showers [138-140], the Fly's Eye detector [141] and its 
successor HiRes [142] established the measurement of the lon- 
gitudinal development of air showers using fluorescence tele- 
scopes and studied the evolution of the shower maximum with 
energy [143, 144]. Currently, two observatories are in operation 
that use the fluorescence technique for the determination of the 
energy scale and for composition studies: The Pierre Auger Ob- 
servatory in the Southern hemisphere [145] and the Telescope 
Array (TA) in the Northern hemisphere [84]. 



The measurement of the longitudinal air shower development 
with fluorescence telescopes relies on the fact that the charged 
secondaries of an air shower excite the nitrogen molecules in 
the atmosphere that in turn emit fluorescence light. Since the 
light yields [146] are proportional to the energy deposited in 
the atmosphere, this observation allows to reconstruct the lon- 
gitudinal development of the air shower as a function of slant 
depth. 

A typical example of a reconstructed energy deposit profile 
of an ultra-high energy air shower is shown in Fig. 9. For this 
particular shower, the full profile was observed and the total 
calorimetric energy could be obtained by simply adding up the 
data points. In general, however, only part of the profile can be 
detected, because the shower either reaches ground or its ris- 
ing edge is obscured by the upper field of view boundary of 
telescope. Therefore, the profile is usually fitted with an appro- 
priate trial function [147] that allows the extrapolation of the 
shower outside of the field of view and to below ground level. 
Popular choices for fitting longitudinal profiles are the Gaisser- 
Hillas function [111] (used e.g. by Auger [148]) or a Gaussian 
in shower age [149] as it was used for the final HiRes analy- 
ses. The calorimetric energy of the shower is then given by the 
integral of the fitted energy deposit profile. 

In addition to the calorimetric energy, the measurement of 
the longitudinal energy deposit profile provides a direct ob- 
servation of the shower maximum. As can be seen in Fig. 9, 
where simulated longitudinal shower profiles are superimposed 
on the measured profile, even on a shower-by-shower basis a 
rough distinction between heavy and light primaries is possi- 
ble by comparing the position of X max . In principle, the full 
distribution of shower maxima for showers with similar energy 
contains the maximum information about composition that can 
be obtained from fluorescence detectors. Given enough statis- 
tics and an exact knowledge of the expected distributions for 
different primaries, it should be possible to extract composition 
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Figure 9: Example of a longitudinal air shower development as measured with 
fluorescence telescopes. Data points are taken from [145] (E = (30 ± 2)EeV) 
and compared to ten simulated [133] air showers for three different primary 
particle types using the hadronic interaction model Epos1.99 [36]. 



groups (see e.g. [150]) similar to what is done for surface de- 
tectors. In the following, however, we will concentrate on the 
first two moments of the X max -distribution, (X max ) and cr(X max ). 

For the determination of the average shower maximum, ex- 
periments bin the recorded events in energy and calculate the 
mean of the measured shower maxima. For this averaging not 
all events are used, but only those that fulfill certain quality 
requirements that vary from experiment to experiment, but all 
analyses accept only profiles for which the shower maximum 
had been observed within the field of view of the experiment. 
Without this condition, one would rely only on the rising or 
falling edge of the profile to determine its maximum, which 
was found to be to unreliable to obtain the precise location of 
the shower maximum. The field of view of fluorescence tele- 
scopes is typically limited to 1-30 degrees in elevation. There- 
fore some slant depths can only be detected with smaller effi- 
ciencies than others, resulting in a distortion of the measured 
X max -distribution due to undersampling in the tails of the distri- 
bution [151, 152]. For instance, a detector located at a height 
corresponding to 800 g/cm 2 vertical depth cannot detect shower 
maxima deeper than 800, 924 and 1600 g/cm 2 for showers with 
zenith angles of 0, 30 and 60 degrees respectively. On top of 
this acceptance bias an additional reconstruction bias may be 
present that can further distort the measured (X max )-values. 

There are two ways to deal with such biases: If one is only 
interested in comparing the data to air shower simulations for 
different primary particles, then the biased data can be simply 
compared to air shower predictions that include the experimen- 
tal distortions. For this purpose the full measurement process 
has to be simulated including the attenuation in the atmosphere, 
detector response and reconstruction to obtain a prediction of 
the observed average shower maximum, (X max ) obs . Another 
possibility is to restrict the data sample to shower geometries 
for which the acceptance bias is small (e.g. by discarding verti- 
cal showers) and to correct the remaining reconstruction effects 
to obtain an unbiased measurement of (X max ) in the atmosphere. 
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Fi gure lOi Elongation rates obtained by a linear fit in lg£ to the X max data 
of HiRes, Yakutsk, TA and Auger above different energy thresholds. Only fit 
results with x 2 /Ndf < 2 are shown. The yellow, solid band is the average 
obtained for HiRes, Yakutsk and TA , the green hatched band indicates the 
average for all four experiments. 



Whereas the former approach maximizes the data statistics, 
the latter allows the direct comparison of published data to air 
shower simulations even for models that were not developed at 
the time of publication. Moreover, only measurements that are 
independent of the detector-specific distortions due to accep- 
tance and reconstruction can be compared directly. 

The HiRes and TA collaborations follow the strategy to pub- 
lish (X max ) obs [130, 132] and to compare it to the detector- 
folded air shower simulations. In the HiRes analysis the cuts 
were optimized to assure an X max -bias that is constant with en- 
ergy, but different for different primaries and hadronic inter- 
action models. The preliminary TA analysis uses only mini- 
mal cuts resulting in energy dependent detection biases. The 
Auger collaboration quotes average shower maxima that are 
without detector distortions within the quoted systematic uncer- 
tainties [153] due to the use of fiducial volume cuts. Yakutsk 
derives X m . dx indirectly using a relation between the slope of 
the Cherenkov-LDF and height of the shower maximum (cf. 
Sec. 3.2). This relation is derived from air shower simula- 
tions and is universal with respect to the primary particle and 
hadronic interaction models [154]. We will therefore assume 
in the following, that the the Yakutsk measurement is bias-free 
and that it can be compared to air shower simulations directly. 

To allow a comparison of the results of these experiments and 
moreover to calculate (In A) using the Epos model (cf. Sec. 3.4) 
which was not used in some of the original publications, we 
correct the (X max ) obs -values of HiRes and TA by shifting them 
by an amount A which we infer from the difference of the pub- 
lished (X max ) obs -values for proton, QGSJetII to the simulated 
values that are obtained without detector distortions: 



(^max)c 



(^max)obs + A 



(27) 
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Figure 11: Left: Shower-to-shower fluctuations of X mllK in the atmosphere as measured by Yakutsk [127] and Auger [131] compared to air shower simulations [133]. 
Right: Width of a Gaussian fit to the truncated X max -distributions as measured by HiRes [130] compared to air shower simulations including detector effects. 



where 

A = <X max (QGS JetII, p)> - <X max (QGSJ ET II, p)) obs . (28) 

The obtained A-values are approximately constant with energy 
for HiRes with (A) = 26g/cm 2 . A decreases with energy for 
TA from A = 24g/cm 2 at 10 183 eV to A = 9 g/cm 2 at 10 19 7 eV. 
Although these values correspond to up to 25% of the predicted 
proton/iron difference, the correction procedure seems justified, 
because both experiments measured a (X max ) obs that is very 
close to the detector-distorted QGSJETlI/proton predictions and 
the A-correction preserves the residuals wrt. this prediction for 

max)corr- 

The (X max )-values measurements from fluorescence detec- 
tors are shown in Fig. 8. At around the region of the ankle 
of the cosmic ray spectrum the measurements are compatible 
within their quoted systematic uncertainties (~ 10 g/cm 2 ) and 
the (X max ) is close to the prediction for air showers initiated 
by a predominantly light composition. Below this energy, the 
few data points from Auger as well as the HiRes/MIA data 
confirm the trend observed by the non-imaging Cherenkov de- 
tectors, namely a large elongation rate indicative of a transi- 
tion from a heavy to a lighter composition. A complete cov- 
erage of the transition region down to 10 17 eV with new fluo- 
rescence detector measurements is expected in the near future 
with the High Elevation Auger Telescopes (HEAT, operating 
since 2009) [155] and the Telescope Array Low Energy Exten- 
sion (TALE, under construction) [156]. 

At ultra-high energies, the experimental situation can be 
quantified by fitting the data with a linear function in the log- 
arithm of energy. The resulting elongation rates are summa- 
rized in Fig. 10. As can be seen, the elongation rates of the 
four experiments agree within uncertainties above a few EeV, 
though the statistical precision deteriorates quickly with in- 
creasing energy. Numerical values for fits above 10 18i eV are 
given Tab. 1. Auger, which has collected by far the largest 
statistics of the four experiments, observes a small elongation 



rate of Dio = 26 + 5 (g/cm 2 )/decade, which could indicate a 
gradual increase of the average mass of cosmic rays at ultra- 
high energies if the elongation rate for a constant composition 
is indeed 50 to 60 (g/cm 2 )/decade as predicted by the mod- 
els. This small value of Dio is confirmed by HiRes, TA, and 
Yakutsk, for which the weighted average of elongation rates 
from Tab. 1 is 27+9 (g/cm 2 )/decade. 

The absolute depths at 10 19 eV are, however, in poor agree- 
ment among the four experiments and the differences of up to 
24 g/cm 2 between the Auger and HiRes results are larger than 
expected for individual systematic uncertainties at the 10 g/cm 2 
level. It is worthwhile noting that without the A correction, the 
three fluorescence detectors agree almost perfectly at ultra-high 
energies. This might be a mere coincidence or a hint that either 
HiRes and TA overestimate their bias or that the assumption of 
A as does not hold for Auger. 

As explained in Sec. 2.1, the fluctuations of the shower max- 
imum provide another composition sensitive observable. The 
measurements of cr(X m . dx ) from Auger and Yakutsk are shown 
in the left panel of Fig. 1 1 . Both data sets were corrected for 
the detector resolution and can thus be directly compared to air 
shower simulations. The Auger data exhibits a significant nar- 
rowing of the X m - dx distributions with energy starting at about 
the energy of the ankle. The low energy width is compatible 
with a light or mixed composition, but the small cr(X max ) at 





X\9 

[g/cm 2 ] 


D l0 

[(g/cm 2 )/decade] 


AT 2 /Ndf 


HiRes 


782+3 


23+11 


1.7/4 


Yakutsk 


773±5 


35+19 


1.9/5 


Auger 


758+1 


26+5 


1.9/5 


TA 


774+5 


32+18 


1.4/4 



Table 1 : Results of a fit with <X max > = X l9 +Di (lg(£/eV) - 19) to <X max >-data 
above 10 186 eV. 
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Figure 12: Comparison of Auger measurements [131] of the average shower depth and its fluctuations with air shower simulations. The energy of the data points 
increases from top to bottom, as indicated by the dashed line. 



high energies points to the presence of a significant fraction of 
CNO or heavier nuclei with little admixture of light nuclei (cf. 
Eq. (12)). The data of the Yakutsk array are compatible with 
a light composition at all energies. Below 10 I93 eV there is a 
good agreement with the Auger results, but their last data point 
is in contradiction with the small width quoted by Auger. 

The fluctuation measurements of HiRes are shown in the 
right panel of Fig. 11. Instead of cr{X m . dx ), HiRes published 
the width of a Gaussian fit to the X max -distributions that were 
truncated at ±2cr(X m - dx ) without correction for detector resolu- 
tion. This variable can then be compared to air shower simula- 
tions including detector effects. As can be seen, HiRes finds a 
large width at low energy that is, similar to Auger and Yakutsk, 
compatible with a light or mixed composition. Above 10 19 eV 
the width remains compatible with proton simulations, albeit 
with large statistical uncertainties that could also accommodate 
a narrower width. 

The compatibility of the (X mdx ) and cr(X m . dx ) measurements 
with air shower simulations can be studied within the cr(X max )- 
(X m . dx ) plane introduced in Sec. 2.1. The Yakutsk data would 
cluster at around the simulated proton values in this plane and 
the HiRes data cannot be analyzed in this way without a full 
detector simulation. Therefore only the Auger data are shown 
in Fig. 12. 

If hadronic interactions at ultra-high energies are modeled 
correctly and if the cosmic ray composition is any mixture of 
elements between proton and iron, then the data points must lie 
within the contours shown in Fig. 12. As can be seen, this is not 
the case for the outdated QGSJetOI model. For both, QGSJetII 
and Sibyll2.1, the Auger data are barely at the edge of the con- 
tour, which would imply a somewhat unnatural transition from 
a proton- to helium- to nitrogen-dominated composition. Using 
Epos1.99, the corresponding composition would be mixed at 
low energies and very nitrogen-rich at high energies. Whereas 
these considerations clearly demonstrate the power of a com- 
bined observation of {X mdx ) and cr(X mdx ), the current system- 
atic uncertainties of the Auger measurements do not allow for 
a stringent test of the models: If the Auger data are shifted si- 
multaneously by ±syst(£'), +sys((X max )) and ±syst(cr(X max )) as 
in [157] (indicated by lines in Fig. 12), then even for QGSJetOI 



a marginal compatibility remains. 

3.4. Average Logarithmic Mass 

As it should have become clear in the previous sections, the 
measurement of the mass composition of cosmic rays from air 
shower data is a challenging task. There are still considerable 
systematic differences between the various air shower measure- 
ments that need to be resolved. A further obstacle arises from 
the large shower-by-shower fluctuations of air shower observ- 
ables which are found to fluctuate at a level which, for light 
primaries, may become comparable to the mean separation of 
proton and iron primaries. Moreover, and maybe most impor- 
tantly, the interpretation of these measurements in terms of the 
chemical composition of cosmic rays relies on the validity of 
models of hadronic interactions at ultra-high energies. 

The best discrimination power in terms of models for the ori- 
gin of ultra high energy cosmic rays is of course given by a 
determination of the energy dependence of different groups of 
primary elements. But as could be seen in e.g. Fig. 5, even 
when different experiments use the same hadronic interaction 
model to interpret their data, there are large differences in the 
quoted fluxes of elemental groups and similar discrepancies can 
be observed when comparing the unfolding results using differ- 
ent interaction models. Part of these difficulties arise due to 
the large negative correlations between the unfolded elemental 
fluxes, which are however not quoted by the experiments. 

We therefore restrict ourselves to a comparison of the average 
logarithmic mass of cosmic rays, which we obtain from either 



(]nA) = Yjfi lnA ' 



(29) 



for experiments that derived elemental groups of mass A, with 
flux fractions /j or from 



(In A) 



In 56 



(30) 



in case of experiments that measured (X mdx ). For the latter, 
model uncertainties can be estimated by using different pre- 
dictions for (X^ ax ) and (X^). This comparison is shown in 
Fig. 13, where we also include the QGSJetOI model for com- 
pleteness, due to the ubiquity of results from particle detectors 
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that used this superseded model. Obviously, the systematic dif- 
ferences in (X m . dx ) discussed in the last section propagate di- 
rectly to (In A). To guide the eye and to be able to compare 
the results from optical detectors with those of particle detec- 
tors (see below), the upper and lower (In A) ranges are sketched 
in Fig. 13 by solid lines. As can be seen, the experimen- 
tal systematics in (X max ) translates to an uncertainty of about 
cr((lnA)) « +0.5. The composition trends that were already 
visible in Fig. 8 can again be observed in (InA): All model 
interpretations suggest a gradual increase of the average loga- 
rithmic mass of cosmic rays between 10 15 eV and 10 17 eV fol- 
lowed by a transition towards a lighter composition during the 
next decade. The heaviest composition with (InA) w 3.5 fol- 
lows from the Tunka data interpreted with QGSJetII at around 
10 17 eV. The (InA) values of HiRes and TA are compatible with 
a pure proton composition when using one of the two QGSJet- 
fiavors. A trend towards a heavier composition would follow 
from Auger data for all models and also for HiRes and TA if 
interpreted using Sibyll or Epos. It is interesting to note that 
the next version of QGSJetII [158] for which some model pa- 
rameters were re-tuned to new data from the LHC will have a 
similar (X mdx ) as Sibyll and thus the combination of any of the 



C^max) data with one of the contemporary versions of the three 
available interaction models will result in a (InA) significantly 
different from zero at ultra-high energies. 

Particle detectors usually do not publish air shower observ- 
ables but directly the interpretation in terms of elementary frac- 
tions, and in that case only the differences between models with 
which the data were analyzed can be used for a limited estimate 
of the theoretical uncertainties. Results that were obtained with 
out-dated interaction models like e.g. the AGASA measure- 
ments [159] will be ignored in the following. Since usually only 
fractions of elemental groups are quoted it is not obvious which 
value of InA, to assign in Eq. (29). To translate the data from 
Tibet ASy [89] into (InA), we assume equal fluxes of protons 
and helium and assign to 'heavy' fragments A = 32. However, 
we note that the chosen procedure of comparing fluxes from 
different measurement campaigns with different event selection 
and energy calibration may introduce additional systematic un- 
certainties particularly in view of the steep power-law spectra 
involved, which we can not account for here. For KASCADE- 
Grande [92], where the intermediate mass group is composed 
of He, C, and Si, we again assume equal fluxes and take the 
logarithmic mean of A 12. For data that were analyzed in 
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a simple bimodal proton/iron model like [99, 101] the (lnA) 
calculation is technically easy, but it is difficult to assess the 
systematic uncertainty arising from this simplified model. Data 
from EAS-TOP are based on electrons and GeV muons [160] 
as measured in the calorimeter at the surface as well as on elec- 
trons and TeV- muons, the latter measured in MACRO [81]. Of 
all the experimental particle detector results studied here, only 
Auger published the measured air shower observables rather 
than their interpretation. Since the average muon production 
depth and the rise time asymmetries are well correlated with 
Xmax we assume that they also depend linearly on (In A) and 
can therefore use the air shower simulations folded with the de- 
tector response from [1 17] and [1 10] to estimate (In A) from the 
equivalent of Eq. (30) for these variables. 

The resulting energy evolution of (In A) as derived from par- 
ticle detector data is displayed in Fig. 14 for different hadronic 
interaction models. The upper and lower experimental bound- 
aries from optical detectors are indicated by the superimposed 
lines. As can be seen, the systematic differences between exper- 
imental results at low energies are considerably larger than in 
the case of optical detectors spanning a range of up to A(ln A) « 
±1. Nevertheless, all experiments below 10 17 eV report a rise of 
(In A) with energy that could be reconciled with the (X m . dx ) re- 
sults by an appropriate rescaling. In the energy region toward 
the ankle, surface detector data are sparse. The Haverah Park 
results tend towards a lighter composition at 10 18 eV, though 
with large statistical uncertainties. At ultra-high energies only 
the surface detector data from Auger are available for an in- 
terpretation with modern hadronic interaction models. For both 
simulations, using QGSJetII and Sibyll2. 1, these data are com- 
patible with an increase of (In A) above 10 19 eV. 



4. Search for Neutral Primaries 

Measurements of neutral primaries, i.e. neutrons, photons, 
and neutrinos provide additional crucial information about the 
acceleration models and sources of cosmic rays as well as on 
their propagation through the universe. Unlike charged cosmic 
rays they are not deflected by magnetic fields and could point 
back to their sources. Specifically, high energy neutrinos point- 
ing back to a source would provide the first unambiguous signal 
of a high energy cosmic ray accelerator. Moreover, a 'guaran- 
teed' flux of high energy photons and neutrinos is provided by 
the GZK-effect (see Sec. 1). An observation of such cosmo- 
genic neutrinos or photons is considered as a 'smoking gun' that 
would complement the observation of the flux suppression seen 
in the cosmic ray energy spectrum. Neutrons arise at all ener- 
gies from charge exchange interactions at the source or in the 
interstellar medium. However, their lifetime limits their propa- 
gation distance to about lOkpc ■ (Zs/EeV). 

Separating neutron- from proton-induced showers is not pos- 
sible. Nevertheless, samples of air showers can be enriched 
with respect to light primaries and searches for neutron point 
sources can be performed. Such studies have yielded no de- 
tection yet and are reviewed in [21]. Separating photon- or 
neutrino-induced EAS from nuclei-induced showers is exper- 
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Figure 14: Average logarithmic mass of cosmic ray as a function of energy 
derived from particle detector measurements for different hadronic interaction 
models. Lines are the upper and lower bounds from optical measurements (cf. 
Fig. 13) for the corresponding model. 



imentally much easier than distinguishing light from heavy nu- 
clear primaries. 

In this chapter we will review recent photon and neutrino 
searches aimed at highlighting their complementarity to mea- 
surements of the nuclear composition. 
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4.1. Photons 

Air showers induced by primary photons develop an almost 
pure electromagnetic cascade and differ from hadron induced 
ones by their low number of muons and their deep X m . dx . This 
is illustrated in Fig. 15 for different primaries at 10 I9 eV. Since 
mostly electromagnetic processes are involved in the shower 
development, the predictions are more reliable and do not suf- 
fer from uncertainties in the hadronic interaction models. The 
Pierre Auger Observatory has the advantage of being a hybrid 
observatory in which the position of the shower maximum is 
measured by fluorescence telescopes with simultaneous access 
to the muon number from the surface detector stations. In case 
muons cannot be identified, such as in pure scintillator arrays 
located at the surface, the particle densities at large distances 
from the shower core combined with measurements of the en- 
ergy by e.g. fluorescence telescopes can still provide good dis- 
crimination power. 

Early searches for UHE photons, summarized e.g. in 
Ref. [169], were motivated by particle physics models of 
UHECR origin. In these so-called top-down models (see 
e.g. [170] for an early review) the highest energy cosmic rays 
are decay products of super-heavy relic particles or topological 
defects (TD) left over from the inflationary epoch and which 
are locally clustered in the galactic halo as cold dark matter. 
Their decay yields a relatively high flux of UHE photons and 
neutrinos reaching up to 90 % of the CR flux. The Z-Burst 
model [171, 172] involves resonant interactions of very ener- 
getic neutrinos with the relic neutrino background. These mod- 
els are compared to experimental data in Fig. 16. The inte- 
gral photon flux limits from AGASA [161], Yakutsk [162], and 
TA [165] are based on muon numbers and curvature of the 
shower front. The hybrid data of the Pierre Auger Observa- 
tory, employing X max and various SD-observables, provide the 
best present limits up to 10 19 eV while the Surface Detector data 
alone extend the best present limits up to 5 x 10 19 eV due to the 
much higher statistics [163, 164] available in this data set. 
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Figure 16: Integral photon flux limits at 95% C.L. from AGASA [161], 
Yakutsk [162], Auger [163, 164] and TA [165] compared to flux predici- 
tons for GZK-photons [166], top-down scenarios of super-heavy dark matter 
(SHDM) [167] and topological defect (TD) models, and Z-bursts [168]. 



The results demonstrate that particle physics motivated top- 
down models are strongly disfavored giving support to an as- 
trophysical origin of UHECR. 

The two bands below the experimental data depict photon 
fluxes expected from interactions of UHECR with photon back- 
ground fields, most prominently the Cosmic Microwave Back- 
ground (CMB) radiation. Both, proton and iron primaries have 
been considered in these calculations. A power law energy 
spectrum is assumed with index y = — 2 and a maximum CR 
energy at the source of 10 21 eV. The fluxes are normalized 
to the number of events measured by Auger [173] for ener- 
gies E > 10 18 4 eV. The upper and lower bounds for each pri- 
mary result from more or less favorable radio background fields 
and different cosmological source evolutions [166]. Accord- 
ing to these simulations, chances of observing GZK-photons 
appear well in reach if the source spectrum at the highest en- 
ergies is dominated by light primaries. Moreover, further im- 
provements, besides increasing statistics of data, may be ex- 
pected from multivariate analyses of fluorescence and surface 
detector data. Also, further upward and downwards modifica- 
tions of the predicted photon flux may result from different en- 
ergy spectra and/or different maximum energies of the nearby 
sources [166, 174]. 

4.2. Neutrinos 

Possible neutrino primaries may be the easiest to identify in 
EAS experiments because of the many orders of magnitude dif- 
ference between the electroweak and hadronic cross sections. 
With a neutrino-nucleon cross section of about ct yN - 10 nb 
at 1 EeV, neutrinos may interact at any point in the atmosphere 
or even in the rock of the Earth. Thus, showers starting very 
deep in the atmosphere are likely to be initiated by a primary 
neutrino. Identification of such showers is easiest in horizontal 
direction (9 > 85°) where the atmospheric thickness is effec- 
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Figure 17: Compilation of 90%CL single-flavor upper limits for diffuse neu- 
trino fluxes assuming a proportion of flavors of 1:1:1 due to neutrino oscilla- 
tions. Data are from IceCube [175], Auger [176] and ANITA [177]. The shaded 
area corresponds to expected GZK neutrino fluxes computed under different as- 
sumptions of source evolution scenarios [166] with power-law energy spectra 
of y = -2.0 and £^ ax = Z ■ 10 20 eV. The grey band depicts different transition 
models and source evolutions adapted from Ref. [178] (see text for details). 



tively 30 times that of the vertical direction. Thus, horizontal 
showers containing still an appreciable electromagnetic com- 
ponent (so-called 'young showers') can only be caused by pri- 
mary neutrinos. Larger experimental apertures than for near- 
horizontal showers induced by so-called 'down-going' neu- 
trinos are reached for 'up-going' tau neutrinos skimming the 
Earth at zenith angles between 90° and 96°. After a possible 
v T -interaction, the produced r will be able to escape the Earth 
and decay in the atmosphere or close to ground producing an 
upward-going EAS. 

Such signatures have been searched for in ground based EAS 
observatories such as Auger [179] and HiRes [180] reaching 
highest sensitivities around 10 18 eV, thus matching very well 
the peak in the differential Eydj Y /dE Y fluxes of cosmogenic 
neutrinos. At even higher energies, radio based observations 
such as ANITA [177] start to take over, while IceCube pro- 
vides unprecedented sensitivities towards lower energies. The 
three experimental approaches complement each other very ef- 
fectively with each of them having improved their sensitivities 
by orders of magnitudes in the last decade (cf. Fig. 17). 

Similarly as for UHE photons, these searches did not detect 
any neutrino fluxes yet but provided upper limits. Fig. 17 de- 
picts some exemplary reference predictions again for sources 
emitting proton and iron primaries [166, 178]. Due to their 
low cross-section, neutrinos — other than photons — suffer 
only adiabatic energy loss, so that they can arrive at Earth 
from cosmological distances. As a result, a large level of un- 
certainty in the flux prediction results from the unknown cos- 
mological evolution of source luminosity, as has been realized 
already in [181]. The upper dashed bands shown in Fig. 17 
are obtained for an FR-II type distribution while the lower one 
results for a star formation like evolution. As for the pho- 



ton predictions, power law source distributions with an injec- 
tion index of y = — 2 have been assumed as a reference with 
Ef n!lx = Z ■ 10 20 eV. A fit to the Auger CR energy spectrum per- 
formed e.g. in [166, 182] yields somewhat steeper source spec- 
tra but higher maximum energies of the sources with essentially 
compensating effects with regard to the neutrino fluxes. The 
grey band, adapted from Ref. [178], depicts neutrino predic- 
tions covering different transition models (dip-model at top and 
ankle model at bottom) and source evolutions. 

Atmospheric neutrinos in the TeV to PeV range are primarily 
of interest, because they constitute the background in neutrino 
telescopes such as IceCube [175] or Antares [183]. The cosmic 
rays producing this 'calibrated' beam of neutrinos most domi- 
nantly originates from energies above the knee. Thus, an inter- 
esting question to ask is the level of uncertainty in the neutrino 
flux calculations that arises from the poorly known composition 
(see e.g. [184]). Expressed differently, one may ask whether a 
comparison of the measured atmospheric neutrino spectra with 
flux calculations based on cosmic ray energy spectra could pro- 
vide independent information about the composition at the knee 
and above. A quantitative study of this question has been pre- 
sented very recently in Ref. [185]. It used the two different 
sets of energy spectra shown in Fig. 5 from KASCADE [65] 
combined with direct measurements of cosmic rays at lower 
energies as input to a full Monte Carlo simulation. Based on 
this study, the authors concluded that uncertainties in the all- 
particle cosmic ray flux appear more important for the neutrino 
fluxes than the actual uncertainties of the composition arising 
from the unfolding based on QGSJet versus Sibyll interaction 
models. 



5. Conclusions 

In 1983, Linsley counted the number of contributions sub- 
mitted to the 18th ICRC that were related to the composition of 
cosmic rays and concluded [51] 

Assuming that people wouldn't spend effort on this 
problem unless they believe it can eventually be 
solved, I think the numbers are encouraging. 

Today, almost 30 years later, the quest for the understanding 
of the primary composition of cosmic rays continues and, al- 
though the problem is not fully solved yet, there has been con- 
siderable progress in at least our qualitative understanding of 
the primary cosmic ray composition thanks to the wealth of 
new data collected by both cosmic ray observatories and par- 
ticle physics experiments. 

At the energies between 10 15 eV and 10 17 eV all experiments 
observe energy dependent changes in the shower development 
that are compatible with increasing average mass of cosmic 
rays. Moreover, unfolded spectra of mass groups in this energy 
range suggest that the change of composition is attributable to 
a consecutive cut-off in the flux of the individual mass compo- 
nents starting with protons at a few 10 15 eV up to iron at around 
10 17 eV (cf. Figs. 5 and 6). These composition measurements 
are thus compatible with an interpretation of the knee in the 



18 



< 




| | data (optical) 

SNR&ankle 
-a- SNR&AGN 
— ±— Hypernovae 
— GRB (cannon ball) 

e + e" dip 

mixed, large E 

a max 

— e— mixed, small E max 
— •— galactic sources 




E[eV] 



Figure 18: Compilation of (lnA)-values from recent astrophysical models compared to the uncertainty range of experimental data (optical detectors, see text). 
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dip model [8], the cannonball model [186], a hypernova model [187], and the SNR-AGN model of [188]. 



particle flux due to a rigidity dependent leakage of cosmic rays 
from the galaxy and/or a rigidity dependent maximum energy of 
galactic sources. At the same time, a particle physics interpreta- 
tion of the knee is disfavored, since the interaction models used 
for the interpretation of cosmic ray data were found to bracket 
particle production measurements from the LHC [62] at 7 TeV 
center of mass energy, corresponding to a primary cosmic ray 
energy of 10 164 eV. 

Towards the energy region of the ankle, air shower measure- 
ments indicate a decrease of the average mass of cosmic rays. 
Both the shower-to-shower fluctuations and the average shower 
maximum are compatible with a predominantly light compo- 
sition at a few 10 18 eV. Estimates of the proton-air cross sec- 
tion [189-191] at this energy agree well with the extrapolations 
used in hadronic interaction models, therefore also here a dras- 
tic change of the interpretation of the measurements in terms 
of cosmic ray composition due to uncertainties in air shower 
simulations seems unlikely. 

At the highest energies, above 10 19 eV, the experimental un- 
certainties are still too large to draw firm conclusions from 
the data. The measurements from Auger of (X max ), cr(X max ), 
muon production depth and rise-time asymmetries may be in- 
terpreted as a transition to a heavier composition that may be 
caused by a Peters-cycle in extra-galactic sources similar to 
what has been observed at around the knee. However, the 
(X max )-measurements from HiRes, TA and Yakutsk indicate a 
systematically lighter composition at these energies but with an 
elongation rate compatible with Auger (cf. Tab. 1). 

These experimental differences together with the uncertain- 
ties of hadronic interactions make an astrophysical interpreta- 
tion of the (InA) estimates at the highest energies difficult. In 
Fig. 18 we present a compilation of astrophysical models for 
the composition of cosmic rays together with the (InA) inferred 



from air shower measurements. Here, we show only results 
from optical measurements, since these are more abundant over 
the full energy range and — judging from the differences of re- 
sults from surface detectors at low energies in Fig. 14(a) — are 
also less affected by experimental systematics. The gray band 
is the maximum and minimum of the envelopes from Fig. 13, 
i.e. the upper and lower level of experimental and model differ- 
ences. The curves represent predictions by recent models about 
the origin of cosmic rays with a focus on models at ultra-high 
energies (see [192] and references therein for a comprehensive 
list of models of cosmic rays around the knee). A fairly good 
description of (InA) over the entire energy range is given by 
the two-component model of [188] in which it is assumed that 
CRs up to 10 17 eV are produced in galactic supernova remnants 
while the dominant component at higher energies is of extra- 
galactic origin produced at the shock created by the expand- 
ing cocoons around active galactic nuclei. The dashed violet 
curve shows the classical ankle model and the two red lines 
the so-called dip-model [8] (actual calculation used are from 
[11]). Other than the ankle model, the dip model assumes that 
the transition from galactic to extragalactic models occurs be- 
low the ankle and that the ankle is caused by e + <?~ -interactions 
of protons in the CMB rather than by the onset of the extra- 
galactic cosmic ray component. To make this Bethe-Heitler 
process work, the composition above 10 18 eV must be domi- 
nated by protons. The blue curve represents the generalized 
cannonball model [186] in which cosmic rays are described as 
being ions of the interstellar medium that encountered cannon- 
balls — highly relativistic bipolar jets of plasmoids originat- 
ing from supernova explosions and GRBs — and were mag- 
netically kicked up to higher energies. In Ref. [187] (shown 
as the full magenta line with triangles) it has been suggested 
that hypernova remnants, with a substantial amount of energy 
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in semi-relativistic ejecta, can accelerate intermediate mass or 
heavy nuclei to ultra-high energies. A heavy composition at the 
highest energies is obtained in a model where UHECRs were 
produced by gamma-ray bursts or rare types of supernova ex- 
plosions that took place in the Milky Way in the past [9]. Fi- 
nally, the models labels 'mixed' with large and small E max , are 
taken from Refs. [11, 23]. Here, a source composition similar to 
galactic CRs is assumed with the ankle marking again the tran- 
sition from galactic to extragalactic CRs. In the first realization 
of the model it is assumed that the maximum energy follows 
£max = Z ■ 10 20 5 eV with an injection index at the source of 
y = -2.0 [11] whereas in the latter £ max = Z ■ 4 x 10 18 eV and 
y = -1.6 [23]. 

While data and models agree reasonably well up to the an- 
kle energy, there is almost the full range of masses covered by 
models at energies above 10 19 eV. This may be understood from 
the large uncertainties of experimental data points in this energy 
range. At present neither a composition dominated by light pri- 
maries nor by heavy primaries can safely be ruled out by the 
data. 

As discussed in this paper, a large part of these uncertain- 
ties can be attributed to the uncertainties of hadronic interaction 
models, but also systematic differences between the experimen- 
tal data itself contribute significantly to the large width of the 
(lnA)-range at the highest energies. Fortunately, hadronic in- 
teraction models will improve considerably in the near future 
due to the new particle production data from LHC. At the same 
time further improvements of the quality of air shower data are 
expected because of better instruments and the prospects of a 
combined analysis of the data of surface and optical detectors 
with either existing hybrid detectors or enhancements such as 
e.g. shielded muon detectors in Auger [193]. 

Additional constraints about the charge of cosmic rays may 
be acquired from the anisotropy (or lack thereof) of the arrival 
directions of cosmic rays. The Pierre Auger Collaboration has 
reported directional correlations of the most energetic particles 
(E > 5.5 x 10 19 eV) with the positions of nearby AGN [194- 
196]. Since the turbulent component of the galactic magnetic 
field randomizes the arrival directions of particles with small 
rigidities (i.e. large charge and therefore mass), such a correla- 
tion seems to suggest a dominance of light particles at these en- 
ergies, but it may be still possible that a large scale anisotropy, 
e.g. from the super-galactic plane, remains anisotropic even af- 
ter the passage of heavy particles through the turbulent mag- 
netic field [197]. 

A next generation giant observatory should be able to collect 
large statistics for both, charged particle astronomy and com- 
position analyses above the GZK threshold. If cosmic rays are 
extragalactic and their flux at earth is not dominated by very 
close sources, then the composition analysis at GZK-energies 
will be facilitated by the suppression of any intermediate nu- 
clei between proton and iron at earth due to the different energy 
loss lengths of nuclei [11]. In such a simple situation with a 
'guaranteed' composition of protons and iron only it should be 
possible to simultanously determine the properties of cosmic 
rays and hadronic interactions at energies above ~ 6xl0 19 eV. 



So without any doubt, measurements of the cosmic ray com- 
position up to 10 20 eV remain one of the most important tasks 
for the future and the key to the understanding of the origin of 
the most energetic particles in the Universe. 
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